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Abstract 

We present a study of the behaviour of metastable vortex states in mesoscopic super- 
conductors. Our analysis relies on the London limit within which it is possible to derive 
closed analytical expressions for the magnetic field and the Gibbs free energy. We consider 
in particular the situation where the vortices are symmetrically distributed along a closed 
ring. There, we obtain expressions for the confining Bean-Livingston barrier and for the mag- 
netization which turns out to be paramagnetic away from thermodynamic equilibrium. At 
low temperature, the barrier is high enough for this regime to be observable. We propose 
also a local description of both thermodynamic and metastable states based on elementary 
topological considerations; we find structural phase transitions of vortex patterns between 
these metastable states and we calculate the corresponding critical fields. 

■ PACS: 74.25.Ha, 74.60.Ec, 74.80.-g 

OO 

oo ■ 1 Introduction 

O , 

A significant amount of work has been recently devoted to the study of aluminum superconduct- 
ing disks in a mesoscopic regime, in which the size of the sample is comparable with both the 
coherence length £ and the London penetration depth A. In a first set of experiments jjj the mag- 
netization of such systems was measured at thermodynamic equilibrium. Beyond the Meissner 
state, a series of discontinuous jumps appear when the applied magnetic field is increased. This 

■ behaviour corresponds to the entrance of individual quantized vortices in the sample. A quantita- 
tive understanding of the phenomena involved can be achieved through a numerical study of the 
Ginzburg-Landau equations ^ ||, ||, ||] . Analytical progress can be made if one considers special 
limits of these equations (such as the linear limit Q , the London regime 0, |8| or the Bogomol'nyi 
dual point In later experiments 10 , it was noticed that when sweeping down the applied 
field the sample exhibits a paramagnetic Meissner effect, i.e. it has a paramagnetic magnetiza- 
tion. Such a behaviour previously found in high T c superconductors |Tl| had been interpreted as 
a special features of these materials. However, its occurrence in aluminum disks calls for a less 
specific explanation which results from the role played by metastable states at temperatures well 
below the critical temperature. 

This work is devoted to the study of metastable states in a mesoscopic superconductor. We 
shall investigate vortex patterns in the London regime, i.e. in the limit c» 1, where k is the 
Ginzburg-Landau parameter. Using known results p"2| for the magnetization of an infinite cylinder 
with a circular cross-section of radius R, we derive, in the mesoscopic regime, closed expressions 
for the magnetic field and the free energy that are suitable for analytical calculations. They allow 
to obtain the Bean-Livingston barrier and to calculate the series of matching fields (critical fields 
corresponding to the entrance of the iV-th vortex) . The paramagnetic Meissner state can then be 
explained using a simple scenario for the hysteretic behaviour of the magnetization, involving the 
appearance and the disappearance of the surface barrier. 
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Once the existence of metastable states established, a criterion to classify them is required. 
The number of vortices is a simple topological number but it is not precise enough to distinguish 
between all different type of states. Projection methods into eigenstates of the linearized the- 
ory produce a set of integers that were proposed as classifying numbers to characterize a given 
metastable vortex state Q and that would generalize quantum numbers to a non-linear theory. 
However, these numbers are non-generic and it is not obvious to relate them to the geometric 
features of a given pattern. A more geometrical path will be followed here to present a topological 
description of vortex patterns. Our study is based on an analogy between a configuration of vor- 
tices and a dynamical system, obtained by interpreting the superconducting current as the local 
(phase-space) velocity of a particle. The Hamiltonian of the associated dynamical system turns 
out to be the magnetic field. We shall show that a configuration of vortices can be characterized 
by the critical points of the magnetic field (or equivalently the points where the current vanishes). 
The number of critical points of a given type (maxima, minima, saddle-points) will provide a nat- 
ural set of topological invariant integers associated with a configuration. This geometrical analysis 
reveals the existence of structural transitions between states with the same number of vortices: 
indeed, the number and the nature of the critical points can vary abruptly when the applied field 
exceeds some specific values, provoking topological phase transitions. For a simple polygonal con- 
figuration of vortices, we calculate the corresponding transition fields. These topological changes 
are in fact best visualized using a 'dual' description of a system of vortices, closely related to 
the paramagnetic Meissner effect: there exists a special closed curve F, first introduced in ||, 
which is everywhere orthogonal to the current lines and which passes through the critical points 
of the system. Therefore any structural change affecting the critical points, will manifest itself in 
a topological change of the shape of Y. This curve thus provides an efficient and elegant tool to 
characterize vortices patterns and to understand qualitatively the topological transitions between 
different configurations. 

The plan of this article goes as follows. In section 2, we recall how the London free energy can 
be obtained as a limit of the Ginzburg-Landau theory. In section 3, we present the solution of 
London equation in the case of an infinite circular cylinder (details of the calculations are given 
in Appendix A and B). Section 4 is devoted to the mesoscopic limit: we calculate the magnetic 
field and the free energy, study the Bean-Livingston barrier and obtain the matching fields. The 
paramagnetic Meissner state is obtained in section 5. Section 6 is devoted to the topological 
investigation of vortex patterns: we show that the critical points characterize a configuration and 
we demonstrate the existence of topological phase transitions. In the second part of section 6, the 
dual description of a configuration of vortices, via the curve F, is introduced to provide a better 
understanding of the transitions between patterns. Section 7 contains some concluding remarks 
and discussions. 



2 The London limit of the Ginzburg-Landau equations 

Consider an infinite superconducting cylinder with a cross-section of area VL in the constant mag- 
netic field h e z along the axis of the cylinder. In the presence of a magnetic field the Ginzburg- 
Landau free energy density per unit length of the cylinder, defined as difference of free energies 
with and without magnetic field T — Fgih) — Fs(Q), reads 




(l-H 2 ) 2 + |(v-zi)v;| 2 )dV (1) 



T = 



where ip — \ip\e lx is the order parameter and h — V x A is the local magnetic induction. The 
integration is performed over the cross-section of the cylinder. A superconductor is characterized 
by two length scales: the London penetration depth A and the coherence length £. The Ginzburg- 
Landau parameter is defined as their ratio k = A/£. In this work lengths are measured in units of 
A\/2 and the magnetic induction in units of the thermodynamic critical field H c — 0o/2v / 2ttA^, 
where 0o = hc/2e is the magnetic flux quantum. The free energy ([[]) is given in units of ^ 2 H^/4ir. 
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The Ginzburg-Landau equations for -0 and h are obtained from a variation of T . They are 
given by 

-(v-iA\ip = 2k 2 V (1 - |V| 2 ) (2) 

Vxh = 2\tp\ 2 (Vx~ A) (3) 

Equation ([}]) is the Maxwell-Ampere equation which could be written as well using the current 
density j = Im yt/}*Vifj — \ip\ 2 A. Outside the superconducting sample, ip = 0. The boundary 
condition on the surface of the superconductor is obtained by requiring that the normal component 
of the current density vanishes (superconductor/insulator boundary condition p3|): 



(V - iAj V = (4) 



here n is the unit vector normal at each point to the surface of the superconductor. 

We study a superconductor in an applied external field h e the relevant thermodynamic potential 
is the Gibbs free energy G obtained from F via a Legendre transformation 

G = F — flh e B (5) 

where the total (dimcnsionless) magnetic induction (in the z direction ) B is given by the averging 
the magnetic induction over the cross-section 

B = ±Jh(r)d 2 r (6) 

where the magnetic induction is always in z direction: h = hz. 

In the normal phase, we have ip = 0, B — h e and Fs(h) — Fjy(h) = ilh 2 /2. Therefore, the 
corresponding Gibbs free energy Gm is given by: 

O/j.2 

G N = F N (h) - Slh e B = F N (0) - (7) 

At thermodynamic equilibrium, the superconductor selects the state of minimal Gibbs free 
energy Q . The quantity which is measured in experiments is the magnetization M of the super- 
conductor due to the applied field given by AttM = B — h e . It is obtained, at thermodynamic 
equilibrium using the thermodynamic relation ||l3f : 

~ M -hw e (8) 

where the difference of the (dimensionless) Gibbs energies Q = G$ — Gjv up to a constant equal 
to the superconducting condensation energy, is given by 

g = t -nh e B + nh 2 e /2 (9) 

The dimensionless ratio n ~ A/£ is the only free parameter to describe the superconducting 
state. It determines, in the limit of an infinite system, whether the sample is a type-I or type-II 
superconductor G|. For n > l/\/2, i-e. for type-II superconductors, it is energetically favorable 
for the system in a large enough magnetic field to sustain normal regions in the bulk which appear 
as vortices. Their magnetic flux quantized in units of </V We shall now study the system of 
vortices in an extreme Type II superconductor, i.e. in the London limit k> 1. 

In this limit, the amplitude of the order parameter \ip\ does not vary and can be taken to one 
except at the position of the vortices where it vanishes. The Ginzburg-Landau free energy reduces 
to 

d 2 r (10) 



^=/(f + l(^- a > 
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Using this equation and \ip\ 2 = 1, except at the position of the vortices, we can write 



V - iAj ip 

so that the free energy is given by the expression 

1 



= V x h 



T 



h 2 + I (V x h] 2 ) ,'rr 



(11) 



(12) 



which in our units is independent of n. To obtain the equation for h, we notice that the phase \ 
is non singular except near a vortex (say at r = 0) where it has the property that 

V X . dl = 2vr (13) 

on any closed contour encircling the singularity at point r = 0. Using Stokes theorem, 

<j>^X-dl = J J dS-V xVx = 2tt, (14) 

we deduce that the z-component of V x V% equals 27r<5(r). The equation for h is then given by 
@ 

V 2 h - 2h = ~ATi8{f) (15) 
This equation can be generalized to the case of N vortices placed at points fk 

N 



V 2 h- 2h = -4vr^ 5{r - r k ) 



(16) 



k=l 



In the next section we shall present the solution of this equation in the particular case of a cylinder 
with circular cross section. 



3 Solution for N vortices in a circular cylinder 

^From now on we shall study a circular cylinder of radius R. The magnetic induction and Gibbs 
energy of a system of N vortices in this geometry has been studied first in |L2]]. The detailed 
derivation of the results is presented in the Appendix A. The solution h(r, 0) of the equation ( |l6| ) 
satisfying the boundary condition h(R, 9) = h e can be written as a sum of three terms: 

h(r, 6) = h M (r) + h v (r, 9) + hy{r, 9) (17) 

The first term 

h M (r) = hJ (V2r)/I (V2R), (18) 

where Iq(x) is the modified Bessel function of first kind, describes the Meissner effect in the absence 
of vortices. The second term is the magnetic induction of N vortices placed at points = (r^, 9^) 
for k = 1, . . . , N and it is given by 

N 

h v {r,d) = 2 y £ j K (V2\?-r k \) (19) 
fc=i 

The third term, written with help of the modified Bessel function of first and second kind /„ (x) , 
K n (x) as 

+ °° — T ( F) \ N — 

hy(r,9) = -2 V K n (V2R) ") V ^ V I n (V2r k ) cosn(9 - 9 k ) (20) 
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ensures the boundary condition h(R, 9) = h e . We are interested in polygonal rings of vortices. For 
this particular configuration r& — d, 8 k — 2i:k/N and the expression for hy(r, 6) can be simplified 

h v (r,6)=-2N £ Jr n (V2fi) /n(v ^^y d) cosNnO (21) 

The corresponding Gibbs energy has also been obtained in jl^] . The expression of the dimen- 
sionless Gibbs energy is given by 

Q = T -h e J h(r) d 2 r + nR 2 h 2 e /2 (22) 

The field hy diverges at the positions of the vortices and the corresponding part of the free energy 
must be regularized. The singular contribution coincides with one calculated by Abrikosov [1S|] 
for the case of infinite system and is given by 

Too = N£ + 2tt Ko (V2|fj - r k \) (23) 

where £ = 2 7r(lnK + 0.081) is the one-vortex energy. 

The regular part of the Gibbs energy is calculated in the Appendix [B|. Using the expressions 
@,@ and @ for the magnetic field, one obtains 

Qn-Qo = N£ + 2^^> (V2|r 3 - r k \) - 2nh e jh (l ] 
j^k k=i\ MV2«)/ 

9 ^ X n (^)/„(^r I )/ Tt ((72r 3 ) 
_ 2?r <L TJ72R) cosn(8 z -8 3 ) (24) 

n=— oo i,j=l n\v J 

where Qq = itRh\ — ^T^Tffj) is the Gibbs energy in the absence of vortices i.e. in the 
Meissner state. For the particular case of vortices distributed on a regular polygon of radius d, 
the previous expression (E4j) simplifies and becomes 



Go 



£ + 2* f> (|2^sin||) - 2nh e (l - 

- 27 riV V ^n(V2fi)/^(^d) 2 (25) 
„ MV2£) 

Similar expressions have been already obtained previously in [ fl2[ || and we shall now discuss 
them in the limit of a small radius R in order to use the particularly elegant structure of the 
corresponding solutions to derive analytical expressions for the field and the energy. 



4 The limit of a small disk 

We shall now consider the limiting case where R <C A for which the corresponding Gibbs free 
energy was obtained in Q. In this limit the previous expressions greatly simplify. Since all the 
distances are given in units of A-\/2, we can expand the expression for the magnetic induction using 
the expansion of the Bessel functions for small argument, so that 

N 

h(z,z) -h e = M\z\ 2 - 1) + 2X! ln 

1=1 



1 - ZjZ 



(26) 
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where <p e = h e R 2 /2 is the flux of the external field through the system and we use the complex 
notations, 

-a » = (*)«"• < 27 > 

together with complex conjugate z and zj.. This expression for the magnetic induction is solution 
of the modified London equation 



& 2 



N 

h(z, z) — h e = — 7r y 8{z — Zi)8{z — Zi) (28) 



dzdz ' ' 

i=l 



which can be obtained from ( |16[ ) for small variation of the field on the scale of R. 

In the particular case, we shall consider from now on, of vortices distributed along a regular 
polygon of radius d, we have — d and for the N vortices sitting on a shell, the angle 9k of 
the k-th vortex is 9k — 2irk/N. Then, defining the real quantity x = d/R < 1, the relation ( pc| ) 
becomes 

1 - x N z N 

h N (z 7 z) - h e = e (|z| 2 - 1) + 2 In — r= ^ 



(29) 

whereas the Gibbs energy (|25|) becomes 
n Ci 

~ =£"-27r</> e (l-.T 2 ) + 27rln(l-a; 2Ar )-27r(iV-l)lnx-27rln7V (30) 



N 

where Go = Tr(j) 2 e /2 and £' ~ 27rln (iJ/f). 

The relation between the radius a; of the polygonal configuration and the external flux cj) e is 
obtained by minimizing Gn{x, h e ) with respect to x at fixed <p e , namely 

JV-1 Nx 2N - 2 , . 

*> = -^- + T^r (31) 
4.1 The Bean-Livingston barrier 

The relation ( |3l| ) has, as a function of the applied flux <p e , either zero, one or two solutions. The 
latter case corresponds for the Gibbs energy (|50|) to the existence of a Bean-Livingston potential 
barrier fll9| for which one solution is stable and gives the equilibrium position of the vortices while 
the second solution is unstable and gives the height of the potential barrier as shown in the Fig. [l]. 
It should be noticed that for more than one vortex (TV > 1), the energy Qn diverges logarithmically 
at the origin x = due to the repulsion between the vortices. 

The potential barrier disappears for low enough applied flux below a characteristic value 4> m i n 
given by the minimum of the function fn(x) in the rhs of ( |3l| ) as represented in Fig. |^. Defining 
y = x 2 , we have 



N-l Ny 1 *- 1 

f N (y) = ^^ + T^r ( 32 ) 

2y 1 — y" 

N-l , N(N-l)y»-* , ivy*-" 
fM = —2jp-+ + (TT^V)2 (33) 



For large values of N, we have y^ in — 1/2N so that 

N-l N-l 



(34) 



2 'V2N 2 

For 4> e < (f> m in, there is no stable solution and then no equilibrium position for N vortices exists. 



G 




Figure 1: Behavior of the Gibbs energy as a function of the position x of the N = 3 vortex shell. 
The maximum which exists for high enough applied field corresponds to the unstable equilibrium 
point for the vortex configuration and gives the height of the Bean-Livingston barrier. 




Figure 2: Graphical solution of equation The curve shows the behavior of </> e (x). The two 

solutions x s and x u correspond respectively to the stable and unstable positions of the vortex shell 
for values of the magnetic field above <p m i n - 
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Figure 3: Dependence of the Gibbs energy as a function of the applied flux <p e . The equilibrium 
behavior of the system is described by the envelop of the curves for GN{h e ). Arrows show the 
direction in which the magnetic flux is varied. 

4.2 Matching fields 

By plugging the expression ( |3l|) into the expression (^), we obtain the behavior of the Gibbs 
energy for N vortices as a function of the applied flux <j) e . This expression is the counterpart, in 
the London limit, of the parabolae obtained at the Bogomoln'yi point || for n = l/y/2. For a 
given value <p e of the applied field, the optimal number of vortices corresponds to the envelop of 
the set of curves GN(4>e) as represented in the Fig. |[ The matching values 4>n are the values 
taken by the flux of the external field at which Gn-i{4>n) = Gn(4>n)- They correspond to the 
transition between configurations with N — 1 and N vortices respectively. 

Using the relation (|3(]), it is straightforward to obtain the field at which the first vortex 
penetrates the system; it is given by 4>i — £'/27r = ln(i?/£). This field is larger than the field 
at which the Bean-Livingston potential barrier disappears. For larger values of N, in the limit 
<t>e ^ fimin, we approximate the relation ( ^l| ) by 4> e — (N — l)/2x 2 so that 

GN - G0 ~N (fc - + (l + m^]-NlnN (35) 



2tt vt T ' 2 V iV-1. 

Equating this to the similar expression for GN-i{4>e) we obtain that the matching fields for N — > oo 
are given by solution of the following equation 



A " J // l + 21n%^ ) -lnJV + l + 0(l/JV) (36) 



2 V N 



We use the following ansatz (f)^ = ^(aN + b\nN + c). Substituting it to ( |36| ) we obtain the 
coefficients 

2a 2a / , a\ 

b= 5, c= l + ^i -- ), (37) 

a — 2 a — 2 V 2/ 

where the coefficient a satisfies the transcendental equation: 

a = l + 21na (38) 

This equation has one obvious solution a = 1 which is compatible with (^). However from the 
exact numerical calculations of the matching fields using the relation (|3l| ) we infer that the second 
root a ~ 3.5129 is realized in the solution for the matching fields. It leads to the formula 

4> N = 1.76N- 2.32 hxN + 2.3201 - 1.76 + 0{{\nN) 2 /N) (39) 
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In the Fig. ^ the values of the matching fields obtained from the computed free energy are compared 
with our prediction. The discrepancy is attributed to the slow convergence due to the term 
0((lnN) 2 /N). 

4.3 Magnetization and paramagnetic Meissner effect 

A paramagnetic behavior for the total magnetization results from the existence of a high enough 
Bean-Livingston barrier and the long life-time of the corresponding metastable states. We propose 
the following scenario for the hysteretic behavior of the magnetization when the magnetic field is 
increased and then swept down. 

If thermal equilibrium is maintained while increasing slowly the external flux <p e , the results 
of the preceding section hold and the magnetization is given by the derivative of the envelop of 
the curves QN{4>e) according to the thermodynamic relation (|J). This behavior corresponds to the 
upper curve in the Fig. |^. 

When the applied field is decreased, the vortices that are already in the system are confined 
by the Bean-Livingston barrier and cannot escape the sample until the flux cj) e reaches a certain 
minimum value at which the barrier can be overcome, for instance, by the effect of thermal 
fluctuations. For the present discussion we use the ultimate criterion of zero barrier as a definition 
of the transition field. 



^From the discussion of section 4.1, we deduced that the potential barrier for N vortices 
disappears for a characteristic field given for large N by the relation (|34|), (j> m in ~ (N — l)/2. The 
number of vortices is then decreased by one at this field and the energy drops as shown by the 
arrows in Fig. ||. The magnetization calculated from (|J) is shown by the lower curve in Fig. ||. It 
shows clearly a paramagnetic behavior, although it disagrees quantitatively with experiments and 
the numerical simulations based on the exact solution of the Ginzburg-Landau equations. This 
discrepancy results from the rather unrealistic London limit we used for systems for which k ~ 0.3 
in order to obtain an analytic expression for the Gibbs energy. We emphasize that the system is 
not at thermal equilibrium but in a metastable state. It is worth noticing that for values of the 
parameter tfii = £' /2tt = lni?/£ ~ 2, the paramagnetic Meissner effect can be observed as well for 
an equilibrium configuration, i.e. for a minimum of the Gibbs free energy. But this corresponds 
to R ~ £ so that the delta-like vortex approximation, crucial for the present discussion, is not 
justified. 
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Figure 5: Magnetization of the disk as a function of flux <j) e for <\>\ — In (R/£) — 3. The hysteretic 
behavior is indicated by the arrows. 

5 Topological phase transitions between metastable vortex 
patterns 

Up to now, we have considered vortex configurations for which the position x of the vortex ring 
was given by ( pl| ) resulting from the minimization of the energy. We shall now relax this condition 
and study the behavior of the metastable vortex patterns as a function of the applied magnetic 
flux <f> e . Physically, this corresponds to a situation where N vortices are pinned by a potential and 
remain on fixed positions. 



5.1 An effective Hamiltonian system 

When varying the field (j> e , transitions in a vortex pattern with fixed number N of vortices will 
appear. To characterize them, we shall study the behavior of the magnetic field h given by (|29|). 
To that purpose, it is interesting first to notice that, in two dimensions, the Maxwell-Ampere 
equation f= V x hz is the Hamilton equation of a system whose generalized coordinates (p, q) are 
the coordinates (x, y) such that j x = x, j v = y and the canonically conjugated momentum: 

{ \ : ■ w 

Then, the magnetic field corresponds to the Hamiltonian of the system and the current lines are 
the phase space trajectories. For a Hamiltonian system, the flow can be characterized by its critical 
(fixed) points at which 

V/i = (41) 

Using (flof), we see that critical point correspond to zero velocity points such that j = 0. 

In such a description, a vortex corresponds to a maximum in phase-space. However there also 
exist in the system other critical points, such as minima and saddle-points. The number Nm = N 
of maxima (i.e. vortices), N m of minima, and N s of saddles are not independent, because of the 
Euler-Poincare relation |22j , namely the Euler-Poincare characteristic x is given by 

X = N M -N s +N m (42) 

The integer x is a topological invariant and is equal to 1 for a disk thus providing a constraint 
between the numbers of the different types of critical points. Hence, for a system with a given 
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Figure 6: Paramagnetic and diamagnetic zones separated by the curve T. 



number of vortices, the difference N s — N m between saddles and minima is fixed. However, each of 
them can vary and the set of these numbers provides a complete description of the topology of the 
vortex configurations, namely each topological phase can be described using this set of integers. 

A different way to encode the position and the distribution of the critical points of the magnetic 
field in the disk is to use a special contour T, introduced in ||. The curve T is defined by the 
condition that at each point f it is perpendicular to the current density j[f ) . The equation of Y 
is then 

*= r ^M (43 ) 
d9 d e h y ' 

To this definition we must add the requirement that Y must pass through the critical points at 
which j[r) — 0. The curve Y has several branches; one of them encircles all the vortices and defines 
a natural boundary between the diamagnetic (Meissner) domain near the boundary and the inner 
paramagnetic domain which includes the vortices (Fig. [|). An example of this curve is shown in 
Fig. []] in the case of a disk with one vortex. 

The study of the critical points and the resulting behavior of Y provides a complete description 
of the vortex states either for thermodynamic or metastable states and of the transitions between 
different sets of topological numbers when varying the applied field. ^From the behavior of the 
curve r as a function of the applied field <f> e , it might be interesting to draw the analogy with the 
pressure exerted on a closed membrane separating two systems. 



5.2 Topological study of a vortex ring 

We now consider the simple case where all the vortices are placed on regular polygon at a distance 
x from the center of the disk. The critical points can be classified (see Appendix C) according to 
their stability properties, obtained from the hessian matrix = didjh. We review the various 
topological structures obtained and the transitions between them. 

For a ring configuration of N vortices located at the distance x from the center the general 
expression ( |26| ) for the field h in circular coordinates r, 9 reads 

1 - 2x N r N cosN9 + x 2N r 2N 
2x N r N cos NO 



Mr 2 - 1) + In ( ^ N K-jnjr—TTTTZm ) ( 44 ) 



At the critical point the partial derivatives of h vanish. The condition d r h = leads to 

N(x ™ 1)r N-i ^cosA^ + l)-(^ + iyv 

(1 - 2x N r N cos N6 + x 2N r 2N ') (x 2N - 2x N r N cos N9 + r 2N ) ( ' 
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while the condition dgh — rewrites as 



Nr N sin N9 ^ a: 2jV )(l r 2N ) 

(l-2a; Jv r Jv cosJV0 + a; 2Jv r 2JV ) (i 2JV - cos iV0 + r 2Ar ) 1 j 

The solutions of equations ( |45| ) and (^) are studied in detail in Appendix C. One main feature is 
that one has to distinguish between the cases N = 1, N = 2 and N > 2 where N is the number 
of vortices. We now describe these three cases separately. 

5.2.1 Metastable configurations for N = 1 

For one single metastable vortex located at distance x from the center there is no other critical 
point for weak flux (f> e < (j) m {Y), with 

<M1) = (47) 
1 + x 

At (f> e = <j> m (l) a critical point appears at the point r = 1, 6 = —n. For a flux slightly above (f> m (l), 
this point splits into a minimum at r < 1, 6 — —ir and two saddle points on the boundary r = 1 
situated symmetrically about the diameter passing through the vortex. As the flux increases the 
minimum moves towards the center and the saddle points approach the point r = 1, 9 = each 
from its side. They coalesce at this point for <p e = <^f(l), where 

Ml) = ^ (48) 
1 — x 

and for larger values of the flux one saddle point appears inside the disk between the vortex and 
the boundary in the interval x < r < 1, 9 = 0. When (f) m (l) < cj> e < 4>m{1)i the definition (E5) 



must be generalized to the case of critical point on the boundary 22 . It becomes 



X = N M -N s + N m - ^N sb (49) 

where N s f, is the number of saddle points on the boundary. Thus the index is again equal one. 

An efficient way to visualize these transitions is examine the changes in the curve T as the 
magnetic flux varies. In the case <p e > 4>m (1), shown in Fig. |7^, the curve T lies inside the disk 
and separates the paramagnetic (internal) and diamagnetic (external) regions. For cf) e = 4>m{1) 
the saddle point reaches the boundary of the disk and for 4> m (l) < <f> e < 4>m{1) the curve V starts 
and ends on the saddle points on the boundary (see Fig. |?]b). Therefore, in this regime, one cannot 
define a paramagnetic region inside the disk and there is no closed curve that encloses a unit of 
flux. There appears an arc of circle on the boundary where the direction of the current flow is 
opposite to the screening current responsible for the Meissner effect. At <f> e = <p m (l) all the disk 
becomes paramagnetic and the curve T disappears. 

5.2.2 Metastable configurations for N = 2 

For low enough field <f> e < cf> m (2), with 



m (2) = 2^-4. (50) 
1 + x z 



the only critical point is a saddle point at the center. For <j> m {2) < <p e < 4>m(2), where 

M2) = 2^. (51) 

there are two minima situated symmetrically on the diameter perpendicular to the line passing 
through the vortices and two pairs of saddle points on the boundary. This configuration is shown 
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Figure 7: The case of one single vortex (N = 1)- Behavior of the curve T for a). 4> e > 4>m(1) and 
b). <f> m (l) <<j) e < <p M (l) 




Figure 8: The case of N = 2 vortices for a). (j> m (2) < (p e < 4> c (2), b). <fr c (2) < 4> e < <?W(2) and c). 
4>e > 0j\/(2). The curve T is shown. The position of the vortices is taken to be x = 0.7 > y2 — \/3, 
therefore <f> c (2) < (pM (2). 

in Fig. ||a. Each pair is symmetric with respect to the line passing through the vortices and moves 
towards the points r = 1, 6 = or 9 = ir. For <fi e > 4>m(2) each pair disappears and gives birth 
to one single saddle point located between a vortex and the boundary at 8 = or tt as shown 
in Fig. |c. 

Two minima inside the disk move towards the center and at the critical value of the flux 
<f>e — 4>c(2), where </> c (2) = (1 — x A )/x 2 they merge with the central point, which becomes a 
minimum for <j) e > 4> c (2). This critical value can be shown to be always greater than (/> m (2), but 
it can be greater or less than 4>m{2) depending on the value of x. For x < \/2 — \fi w 0.52, one 
has 4> c {2) > 4>m(2), while for x > \/2 — y/3, we have (f> c {2) < 0m(2). The latter case is shown in 
Fig. |b and 

5.2.3 The case N > 2: topological phase transitions 

In that case, the central point is a minimum for all values of 4> e . In addition for low enough external 
flux there are N saddle points inside the disk situated symmetrically at 9 = 2ir(k + 1/2) /N. As 
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Figure 9: The case of N = 3 vortices for a). <f) m (3) < 4> e < </>m(3), b). </>a/(3) < <j) e < <^c(3) and 
c). (f> e > <t>c (3). The curve T is shown. 



the external flux is increased through <f) m (N) , where 

1 - x N 

<P m (N)=N IT - w . (52) 

N pairs of saddle points on the boundary appear together with TV minima in the bulk (see Fig. |^a). 
The minima move towards the center along the line defined by 8 = 2ir(k + 1/2) /N, k = 1, . . . , N, 
while two saddle points constituting the fc-th pair move symmetrically to the points 8 — 2nk/N and 
8 = 2n(k + 1)/N as the flux increases in the interval (f> m (N) < <f) e < 4>m(N). When <fr e = </)m(N), 
with 

l + x N 

4>m(N) = Nj—^ (53) 

pairs of saddle points meet at 8 = 2irk/N and for larger values of the flux there are N saddle 
points situated between a vortex and the boundary. In this regime there are N + 1 minima and 
2N saddle points inside the disk as shown in Fig. [|b. 

At (f> e = 4> C {N), where <ft c {N) is given by the formulae ( |76| ) and (|77|), each of the N minima 
coalesce with each of the N saddle points at r — p X J N ' , 8 = 27r(fc + l/2)/N, where p c is given by 
((76|). For (f> e > (f>c{N) there remains only one minimum in the center and N saddle points. It can 
be shown that <p c (N) is always greater than 4>m(N) for N > 2. 

Then, for <j) e = cj> c (N) a topological transition occurs, when a minimum and a saddle point 
coalesce. This transition and the corresponding curve T is depicted in Fig. ||b and [|c. It is seen 
that the curve V undergoes itself a topological transition. For </> e < 4> C (N) it does not touch the 
center and has a topology of a ring, while for cf> e > (f> c (N) it has a form of petals of a flower 
emanating from the center of the disk. 



5.3 A dynamical interpretation of the curve T 

The curve T, defined above as being everywhere normal to the current flow, has a simple interpre- 
tation in terms of an effective dynamical system, which gives an alternative method to compute 
it. Instead of the hamiltonian flow ([Io|), consider the lines of V/i, obtained from those of J* by 
a rotation of 90°. This flow is clearly non-Hamiltonian and possess fixed points which are not 
present in the hamiltonian dynamics, e.g. the positions of the vortices appear in this picture as 
a point-like sources, though a saddle remain a saddle in this transformation. Then, the curve T 
represents a limit cycle of the flow V/i emanating from the saddle point. This is illustrated in 
Fig. [lO] where it is clearly seen that the same curve as in Fig. |^ a) appears as a limit of the lines 
of V7i. 
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Figure 10: The curve V as the limit cycle of the flow V/i. The parameters are the same as in 
Fig. |a. 



6 Conclusion 

In this paper we have presented two different aspects of vortex patterns in a mesoscopic super- 
conductor. First we focused on equilibrium properties of a type 11 superconductor in the London 
regime. In this case, the Ginzburg-Landau equations reduce to a single linear elliptic equation 
that can explicitly be solved in simple geometries. We have derived closed expressions for the 
free energy and the magnetization of a mesoscopic system which have been obtained as limiting 



cases of a general calculation valid for a system of any size 12 and not from specific assumptions 
related to the mesoscopic regime With these formulae, we have been able to study precisely 
the matching fields and to compute the paramagnetic Meissner effect, resulting from the existence 
of metastable states. The second part of the paper has been devoted to the general and topologi- 
cal study of the metastable states using a parallel with the theory of two dimensional dynamical 
systems. We showed that a vortex pattern can be characterized by as set of topological integers, 
the number of vortices being one of them. Any change in this set corresponds to a topological 
transition of the vortex pattern. Wc calculated the critical fields associated with these transitions 
and we gave a dual interpretation of these transitions. A key concept is the introduction of a spe- 
cial curve r, that embodies the main geometric features of a configuration of vortices. This curve, 
which appears mathematically as a limit cycle of the system of currents generated by the vortex 
pattern, has a simple physical meaning: it separates paramagnetic from diamagnetic domains. 

We emphasize that this approach that focuses on topological aspects are generic and can be 
extended, at least in a qualitative manner, to systems beyond the London regime and with different 
shapes. Here, we considered a circular cylinder to obtain simple analytical expressions. In recent 
studies, analytic expressions for the free energy of a cylinder with an elliptic cross-section were 
found in the mesoscopic limit p8| and, more generally, an electrostatic analogy was developed that 
enables to study mesoscopic superconductors in the London regime with the help of conformal 
transformations [ p9[ ; using such tools, it would be interesting to carry on a quantitative study of 
topological phase transitions in arbitrary domains. 

In this work, only static configurations of vortices have been investigated. A natural extension 
of the concepts introduced here would be to consider the dynamics of a vortex pattern, where a 
vortex would move in the field generated by the other vortices. In this context, the study of the 
deformations of the curve T would certainly shed some light on the mechanism of vortex nucleation 
in a superconducting system. 
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A Magnetic induction in a cylinder 



We present the solution of the London equation ( |15[ ) for the magnetic induction 

N 

V 2 h(r)-2h(r) = -4tt^ £(f - f fe ) (54) 

k=l 

with the boundary condition h(R, 9) — h e . The polar coordinates are convenient for this geometry. 
The equation ( |54| ) with a vanishing right hand side describes the Meissner regime in the absence of 
vortices. The corresponding solution of the homogeneous equation is given by the ^-independent 
function _ 

^>-^ < 55 » 

where the functions I n (x), K n (x) are modified Bessel functions of first and second kind of order 
n which provide the regular and singular solutions respectively of flsi]) . The solution of the 
inhomogeneous equation (Q) for N vortices is written as a sum of the solution of the homogeneous 
equation and a particular solution: 

h(r, 9) = h M (r) + h v (r, 9) + h v (r, 9) (56) 

where hy (r, 9) is the solution for N vortices in the unbound domain 

TV 

h v (r,e)=2^2K (V2\f-f k \) (57) 
fe=i 

and the term hy(r, 9) is introduced to take care of the boundary conditions: 

N 



k=l 



h v {R, 9) = -h v (R, 6) = -2j2 K o V 2 ( R2 + r fc _ 2Rrk cos ^ _ ( 58 ) 



On the other hand, the field hy can be written as a superposition of solutions of (|5j) regular in 
the interior of the boundary: 



hy(R,6)= J2 c n I n {V2r)e me (59) 

n— — oo 

where, since the field is real, the relation c_ n = c* holds. Expanding the lhs of ( |58| ) into Fourier 
series using standard identities of Bessel functions (2?[] we get 



Ko[j2{R2+rl-2Rr k cos(d-O k )})= £ K n (V2R)I n (V2r k )e in ^- e "\ (60) 



71— — OO 



Equating the coefficients of the expansion of both sides of equation (|58| ) one obtains 

= K n (V2R)I(V2r k ) 9k 

I n (V2R) V 7 

from which the expression (|2tj|) follows immediately. 
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B Expression of the Gibbs energy for a cylinder 

One takes advantage of the vector identity (V x h) ■ (V x h) = V • (h x V x ft) + /i • (V x V x ft) 
to rewrite (V X ft) ■ (V X ft) = V • (ft X V X ft) — hV 2 h, since V- ft = 0. Then one uses the London 
equation ([[5]) to rewrite the expression ( [l2|) for the free energy as 

N 

J" = ttV / d 2 rh(f)S(f-r k ) + - / d 2 rW ■ (hxW x h) (62) 



fc=i 



The singular part is obtained from the contribution of hy to the first integral. The regular 
contribution to the free energy can be written using Stokes theorem: 

^-^oo=7rVMr fc )+7rVMr fe ) + - / d6 h(R,d)h r (R,9) (63) 
fc =i fc=i 4 • /o 

where /i r denotes the partial derivative of /i(r, 0) with respect to r and the fields ftjw and h v are 
given by ((lj) and (|2p|). In the last term h{R, 9) — h e can be taken out of the integration sign and 
the integral becomes 

S T n (i/9.r,.~\ / ~ ~ „ _ s 

(64) 



dR OR J 



= 2^ft e |^|| - 4^2^ (lai^K^R) + K Q (V2R)h(V2R)) 

To prove the last equality the properties of the Bessel functions Iq(x) — I\{x) 1 K' (x) = —K 1 (x) 
and the expansion (^) were used. Only n = term contributes to the angular integral. Using 
another property of Bcsscl functions Iq{x)K\{x) + Kq{x)I\{x) = 1/x and one obtains that: 

±* r d e hr (R,9)= ^IhWm ^y^M^) (65) 



We notice that the contribution from the Meissner field Km in (|62j) exactly cancels with the second 
term in the last equation. Hence, the regular part of free energy becomes: 

The term containing the average magnetic induction over the cross section of the cylinder is 
evaluated along the same lines, using the London equation ([l5|): 

k I msr - + 1 1 v . vwy, . v^im + £ (1 - f^) 

(67) 

Using the expressions d66| ) and (67) in the definition ( |22] ) of the Gibbs energy completes the 
derivation of (1241) • 



C Study of the critical points 

In order to study critical points of the Hamiltonian system defined in section 5.1, one has to solve 
the equations (E5|) and (|46|) for a given number N of vortices. 
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Figure 11: Graphical solution of the equation ( |69| ) for one single vortex (N = 1). The position of 
the vortex is taken to be x = 0.5. 

C.l Existence and location of the critical points 

We first look for critical points which are neither on the boundary of the disk nor at the center, 
i.e. we concentrate on the generic position < r < 1. When < r < 1, one deduces that the 
equation j4q ) is satisfied if sin NO = 0, i.e. when = kn/N for k = 0, 1, . . . , 2N — 1. These values 
of the angle can be subdivided into two groups according to the sign of cos NO. We now discuss 
the solution of the equation ( |45| ) for the two cases separately. 

C.l.l Bulk critical points such that cos NO = +1. 

In this case, we are looking for 'parallel' critical points that lie in the same direction as the vortices. 
And we shall prove that for any N, there exists a critical field (/)m(N), defined as follows, 

4> M (N) = N(l + x N )/(l-x N ) (68) 

such that when <p e > 0m (N) there is always a critical point between each vortex and the boundary, 
and no such point exists if <p e < 4>m(N). 

Using the new variables £ = x N and p — r N the equation ( [l5| ) reduces to 

*lp*-i = ^— - (69) 

Looking at this equation, it is clear that systems with iV=l,iV = 2orA r >2 should be discussed 
separately. 

• For N = 1, the plots of the left and right hand sides of this equation are represented in 
Fig. [ll]. They cross for p < 1 only if the rhs evaluated at r = 1 is less than <f> e , i.e for 
e > 4>m = (1 + 0/(1 ~ = (1+ x)/(l — x). In this case the solution of ( |69[ ) can be shown 
to exist using the continuity of the functions. The position of the critical point corresponding 
to this solution satisfies £ < r < 1 , i. e the point resides between the vortex and the boundary. 
For 4> e < (j)M there is no solution. 

• 



For N > 2, it is worthwhile to consider the reciprocal of (|69j) which rewrites as 

^p 1 -^ = -&0(£-p) (70) 
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Figure 12: Graphical solution of the equation ([70]) for N > 2 vortices. The position of the ring of 
vortices is taken to be x = 0.5. 

The right hand side is a convex parabola with zeros at £ < 1 and 1 /£ > 1 . For N — 2 the lhs 
is a constant function and it is clear that if <fi e > 4>m = 2(1 + £)/(l — £) = 2(1 + .t 2 )/(1 — a; 2 ) 
there is always a solution for £ < p < 1, therefore there exists a critical point between 
the vortex and the boundary, while <j) e < 4>m no solution exists. For N > 2 the lhs is 
an increasing fractional power of p (see Fig. |l2|). The equation ([70]) always has a solution 
between £ and 1 if <f) e > cj> M = N(l + £)/(l - £) = ^(1 + x Ar )/(l - a^) and no solution if 
4>e < 4>M- 



C.1.2 Bulk critical points such that cos N6 = — 1. 

We are now looking for 'anti-parallel' critical points that lie in a direction bisecting the angle 
between two neighboring vortices. Here the discussion is more involved: 

• when N = 1 there appears a critical field <fi m (l) — (1— x)/(l+x) such that: when 4> e > </> m (l) 
there exists a critical point in the direction opposite to the vortex whereas when <p e < ^> m (l) 
no such point exists. 

• when N = 2, there are two critical fields </> m (2) and ^> c (2) such that: when m (2) < <j) e < 
</> c (2) there are exactly two anti-parallel critical points and there are no such points when 

4>e < (2) Or (j) e > (j) c {2). 

• when N > 2, one has again two critical fields 4> m (N) and <fi c (N) such that: when (f> e < (j> m (N) 
there exist N anti-parallel critical points, when <fi m (N) < <f> e < <fi c (N), there are 2N such 
points, and when <j) e > 4> C (N), there are no anti-parallel critical points. 

We now prove these assertions and calculate the corresponding critical fields. 

Consider first that N = 1; then the equation (45) in terms of the variables p and £ becomes 

^ = rafci (71) 

And a solution exists in the interval < p < 1 if the value of the lhs at p — 1 exceeds the value of 
the rhs at the same point, as shown in Fig. |l^. This happens for (f> e > m (l) = (1 — x)/(l + x); 
Now, when N > 2, the equation (^5|) can be rewritten as 

^p 1 -^ = r ^(l + £p)(C + p) (72) 
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Figure 13: Graphical solution of the equation ( |7l| ) for one single vortex (N = 1). The position of 
the vortex is taken to be x = 0.5. 




0.5 
P 



Figure 14: Graphical solution of the equation ( |72|) for N = 2 vortices. The position of each vortex 
is taken to be x = 0.5. 



• TV = 2: this case is illustrated in Fig. |J. For (j> e < 4> m (2) = 2(1 -£)/(l +£) = 2(1 - x 2 )/(l + 
x 2 ) there is no solution, for <j) e = <j) m (2) one solution appears at p = 1 and it persists for 



the range of the flux <j) m (2) < 4> e < <^ c (2) 
4> e > ^ c (2), this solution disappears at p — 



2(1 - = 2(1 - x 4 )/x 4 . For higher flux, 



0. 



• The case N > 2 is the most interesting: Fig. 15 shows the behavior of both sides of the 
equation ((72|). Comparing the value of the lhs at p = 1 with that of the rhs and using the 
fact that has a negative curvature, while that of the rhs is always positive, we can 

state that for </> e < <j) m (N), where 



<M7V) = N(l - 0/(1 + = N(l - *")/(! + x iV ) 



(73) 



there is only one solution. Another solution appears at p = 1 for e = ip m (N) and these two 
solution move towards each other when the flux is increased in the interval 4> m (N) < (f> e < 
4> C (N), whereas at the critical value of the flux tp e — <fi c (N) these two solutions coalesce. For 
4>e > 4>c{N), the equation ([72) has no solutions in the interval < p < 1. 
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2|) for N > 2 vortices. The position of each vortex 



Figure 15: Graphical solution of the equation 
is taken to be x = 0.5. 

Calculation of the critical flux <fi c : 

For <j> e (N) = c/> c the function f lhs (p) = and f rhs {p) = (1 + + p)/(l - £ 2 ) are 

tangent at some point p — p c as shown in Fig. [l5|. Equating the value of these functions and the 
value of their derivatives at p = p c we obtain the following equations: 



f upc+m+tpc) = N(i-e)pl ~ 2 
\ M^pc+e + i) = (N-2)(i-e)Pc* 

Dividing the first equation by the second leads to the quadratic equation for p c : 

(N + 2)Zp 2 c + 2(e + l) Pc - (N - 2)£ = 
which has the unique positive solution 

The critical flux is obtained by substituting the value of p c into the first equation of (|74|): 

_ jv(i-e 2 )^"^ 

(Pc + 0(l + ^c) 



(74) 



(75) 



(76) 



(77) 



We notice that for N — 2 the formula (76) yields the correct value p c = 0. Although topologically 
the cases N — 2 and N > 2 are different, the critical flux </> c (2) = (1 — £ 2 )/£j as we shall see later, 
comes out correctly. 

C.1.3 Central and boundary critical points 

So far we did not discuss the two cases, r = and r = 1. For iV = 1 the center, r = 0, of the disk 
does not satisfy the equations (^B|),(46) and therefore is not a critical point. For N > 1 the value 
r = is always a solution of the equations (|45|) , (^6|) . For the points on the boundary, r — 1, the 
equation (^6|) is always satisfied and substituting the value r = 1 into ( ^5|) we find that there are 
2iV critical points on the boundary with an angle given by the following relation: 



cos N9 



2x N 



1 



„2N 



N(l 



JiVn 



(78) 



which is meaningful only for values of the flux between <j) m {N) and 4>m(N). 
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C.2 The nature of the critical points 



We now discuss the topological nature of the critical points and classify them by linearizing the 
stability exponents near each of these points. The problem is then reduced to the study of the 
Hessian matrix of the second derivatives 

<79) 

Since at the critical points inside the disk we have sin NO = the mixed derivative vanish and the 
Hessian matrix is diagonal in polar coordinates. The second derivative with respect to is given 
at the critical points by 



for cos NO = +1 and by 

(i~ P 2 )(i-a 

(i + £p) 2 (£ + p) 2 

for cos NO = — 1. The value of the second derivative with respect to r can be obtained by looking 
at the way the first derivative d r h changes sign near each solution of the equation (^5|). We shall 
now discuss separately the nature of the critical points in the 'parallel', 'antiparallel', 'central' and 
'boundary' cases. 



C.2.1 Bulk critical points with cos NO = +1. 

For (fi e > ipM, there is only one critical point at p\ in the interval £ < r < 1. The results for the 
second derivatives at this point are summarized in the following table: 





< p < pi 


Pi < P < 1 


d r h 




+ 


d' 2 r h 


- 


- 


9jh 




Crit. point 


Saddle 



C.2. 2 Bulk critical points with cos7V6> = -1. 

The following tables provide a classification of critical points in a variety of regimes. 

N = 1 (f> e < (fi m No critical points. 

<fi e > (fi m One critical point at p\ 





< p < pi 


Pi < P < 1 


d r h 




+ 


d' 2 r h 






&ih 


+ 


Crit. point 


Minimum 



N = 2 (fie < (fim No critical points. 

(fim < (fie < <fic One critical point at pi 





< p < pi 


pi < p < 1 


d r h 




+ 


d?h 


- 


- 


Sjh 


+ 


Crit. point 


Minimum 
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<f> e > <p c No critical points. 
N > 2 <p e < <j) rn One critical point at p\ 





< p < pi 


Pi < P < 1 


d r h 


+ 




d?h 




& e h 


+ 


Crit. point 


Saddle 



< (Pe < 



Two critical points at p\ < pi 





< p < pi 


p\< p< pi 


Pi < P < 1 


d r h 


+ 




+ 






+ 


d$h 


+ 


+ 


Crit. point 


Saddle 


Minimum 



> 



No critical points. 



C.2.3 Central and boundary critical points. 

Consider now the point r = 0. Expanding the field h{r, 9) given by (flj) for small r we obtain 



h(r,0) w -{(p e - 2ATlnx) + e r 2 + 2cos(N0) 



i_e 



(82) 



As we mentioned before, for N = 1 the point r = is not a critical point (due to the presence of 
a term linear in r ). For N — 2, two last terms are quadratic in r and the character of the critical 
point depends on the value of the flux <p e . For <p e < (f> c = (1 — £ 2 )/£ the center, r = 0, is a saddle, 
while for <p e > <p c it turns to be to a minimum. This is to be compared with the prediction of 
formula [77]. It is interesting to notice that this happens when the minimum located along the line 
cos 28 = — 1 disappears at the center. For N > 2, the last term in the equation (82) is negligible 
and the center is always a minimum. 

When <j) m < <p e < (pM we must consider the critical points on the boundary. Since for r = 1 
the condition dgh — is satisfied for all 9, we obtain dgh = and d r dgh ^ if cos 2 NO ^ 1, which 
is always satisfied for (p e in the interval between (p m and (pu- Thus the determinant of the Hessian 
matrix (|79|), given by detH — —(d r deh) 2 , is negative and the eigenvalues are of opposite sign. We 
conclude that the critical points at the boundary are saddles. 
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